Supersolid phase and magnetization plateaus observed in anisotropic spin— 3/2 

Heisenberg model on bipartite lattices 
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We study the spin-3/2 Heisenberg model including easy-plane and exchange anisotropies in one 
and two dimensions. In the Ising limit, when the off-diagonal exchange interaction J is zero, the 
phase diagram in magnetic field is characterized by magnetization plateaus that are either transla- 
tionally invariant or have a two-sublattice order, with phase boundaries that are macroscopically 
degenerate. Using a site factorized variational wave function and perturbational expansion around 
the Ising limit, we find that superfluid and supersolid phases emerge between the plateaus for small 
finite values of J. The variational approach is complemented by a Density Matrix Renormalization 
Group study of a one-dimensional chain and exact diagonalization calculations on small clusters of a 
square lattice. The studied model may serve as a minimal model for the layered Ba2CoGe207 mate- 
rial compound, and we believe that the vicinity of the uniform 1/3 plateau in the model parameter 
space can be observed as an anomaly in the measured magnetization curve. 

PACS numbers: 75.10.Kt, 75.30.Gw 67.80.kb 



I. INTRODUCTION AND THE MODEL 

Finding systems - both theoretically and experimen- 
tally - that exhibit novel quantum phases, amongst them 
supersolid states, played an important role in the study of 
strongly correlated systems in the last fifty years. Super- 
fluid (as well as superconducting) phases and quantum 
crystals can be characterized by off-diagonal long range 
order (ODLRO) 1 and diagonal long range oder (DLRO) 
respectively. This classification allows us to think about 
supersolid phases as states in which ODLRO and DLRO 
coexists. Supersolid phases were first observed in the 
context of strongly interacting bosons of 4 He that can 
simultaneously Bose condensed and order in crystalline 
solid^r— Experimental evidence^— of the existence of 
such phase was found after almost half a century, re- 
viving the theoretical interest in supersolid states, and 
indicating that theoretical interpretation might be more 
difficult than the first ideas £r— 

Apparently various bosonic lattice models provide a 
better understanding of supersolid phases. Quantum 
Monte Carlo (QMC) simulations for hard-core bosonic 
Hubbard model on square lattice with nearest-neighbor 
and ncxt-ncarest-neighbor interactions suggested that 
the 'checkerboard' supersolid phase is thermodynami- 
cally unstable, however - through continuous phase tran- 
sition from the superfluid state - a stable 'striped' su- 
persolid emerges^ Similar QMC simulations of a soft- 
core boson model of square lattice indicated a supersolid 
phase that is stable against phase separationj 14 ' 15 

Matsuda and Tsuneto, and Liu and Fisher showed that 
the bosonic picture of supersolid state can be mapped 
onto a model of magnetic supersolid where the magnetic 
order breaks spin rotational symmetry and translational 
invariance at the same timei 16 i 17 Such magnetic analogs 
of supersolid state were observed in triangular lattice via 



Quantum Monte Carlo (QMC) simulations 1 ^— where 
frustration and order-by-disorder mechanism plays an 
important role in the emergence of supersolid phase. 
Classical Monte Carlo simulation on triangular lattice 
supported by mean-field calculation and Landau theory 
suggested that strong anisotropy can stabilize supersolid 
phases^ 1 " Amongst quasi two dimensional systems QMC 
simulations on bilayer dimer models 2 ^— and orthogo- 
nal dimer models^ were also found to exhibit supersolid 
states that are stabilized by strong frustration and/or 
anisotropy. Supersolid states have also been reported in 
the spin-1 Heisenberg chain with strong exchange and 
uniaxial single-ion anisotropies £2r£2. Furthermore a su- 
persolid phase was found in three dimensional spin and 
hard-core Bose-Hubbard model as well.— 

In this paper we investigate spin-3/2 (quantum) an- 
tiferromagnetic models on a square lattice and on a 
chain with both easy-plane and exchange anisotropies 
described by the following Hamiltonian: 



(i,3) (i,3) 

+Aj2(st) 2 + hJ2S? (1) 



where indicates nearest neighbor sites. Our model 
is inspired by the quasi two dimensional Ba2CoGe207, 
where the magnetic spin-3/2 Co 2+ ions form layers of 
strongly anisotropic square lattices 

The paper is structured as follows: In Sec.|H]we discuss 
the phase diagram in the Ising limit and the instabilities 
of the plateaus using perturbation theory. In the fol- 
lowing section (Sec. we map out the phase diagram 
using a variational approximation in different cases, and 
determine the stability of the plateaus and of the super- 
solid phases. To check the reliability of the variational 
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method, we calculate the phase diagram for the spin-1 
model and compare it to the known results in the litera- 
ture. In Sec. HVI a one dimensional chain is studied using 
a variant of the Density Matrix Renormalization Group 
method and evidence for the existence of an intermediate 
supersolid phase is presented. In Sec. [V]we show results 
of an exact diagonalization study on a square lattice. We 
conclude with Sec. [VlJ 



II. THE ISING LIMIT AND 
PERTURB ATIONAL EXPANSIONS AROUND IT 

A. The Ising limit 

The existence of the gapped phase in our model is 
due to the anisotropic terms, so turning off the SfSJ + 

SfSj off-diagonal Heisenberg term what remains are the 
plateaus. For brevity, we call this J —> limit the Ising 
limit. Since the lattice is bipartite and we have nearest 
neighbor interactions only, the spins are not frustrated 
and all the ground states in the plateaus are either uni- 
form or two-sublattice ordered. The ground state wave 
functions and their properties are listed in Table U and 
the phase diagram as a function of magnetic field and 
single- ion anisotropy is outlined in Figure [T] 

Two uniform phases appear in finite magnetic field: 
the fully saturated state with S z = +3/2 on each site and 
the m/m sa t = 1/3 plateau state with S z = +1/2 on the 
sites. We denote these states as P3 and PI, respectively. 
The two-sublattice states include the two antiferromag- 
netic Ising-like states A3 and Al with staggered magne- 
tization \S^ — Sg\ —3 and \S\ — S Z B \ = 1 and vanishing 
uniform magnetization. In addition we find two other 
plateaus, PI and P2, with magnetization that is 1/3 and 
2/3 of the saturation magnetization, respectively. 

The phase boundaries between different phases are es- 
tablished by comparing the ground state energies. A first 
order phase transition occurs between the Al and A3 
phases at A = £J z /2, when the lowest lying energy lev- 
els cross. (£ stands for the coordination number.) The 
ground state degeneracy (4) at the phase boundary is 
just the sum of the degeneracy of the phases it separates 
(2+2). Since the other states are separated by a gap, we 
expect that the level crossing will persist even for finite 
values of J. The phase transition between the phases PI 
and PI is of similar kind. 

The phase boundaries between two-sublattice A3 and 
PI states is more interesting: the ground state at the 
phase boundary is macroscopically degenerate, and goes 
as 2 x 2 JV / 2 . This degeneracy is understood in the fol- 
lowing way: as we cross the boundary by increasing the 
field, the S z = +3/2 spins on the B sublattice do not 
change, while the S z = —3/2 spins become S z = —1/2 
on the A sublattice. At the boundary, the energy of hav- 
ing an —3/2 or —1/2 is equal, thus they create the 2 N / 2 
fold degenerate manifold (N/2 is the number of sublat- 
tice sites). The additional factor of 2 comes from the 



TABLE I: (color online) Summary of ground states in the 
Ising limit. The relevant order parameters in the Ising limit 
are the magnetization m z = \(Sa + S%) and the staggered 
magnetization mj = \\S\ — Sg\. We denote the fully and 
partially polarized antiferromagnetic states by A3 and Al, the 
fully and partially polarized ferromagnetic phases by P3 and 
PI, and finally the plateau states by P2 and PI correspond- 
ing to the 2/3 and 1/3 plateaus respectively. Although, the 
partially polarized ferromagnetic state PI is a plateau with 
m/m aa ,t = 1/3, we (prefer to) call it ferromagnetic state and 
refer to the plateaus as states that exhibit both finite m z and 
m% . C is the coordination number of the (bipartite) lattice. 
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choice of the sublattice (A or B). Turning on J, this 
degeneracy will immediately be lifted (we may think of 
a pseudospin-1/2 Heisenberg like effective model to de- 
scribe this problem), and a gapless phase appears. The 
same scenario holds for the phase boundary between the 
phases PI and P2. These phase boundaries are shown 
by thick red line in Fig. [T] 

Lastly, we examine the phase boundary between the 
uniform and two-sublattice states. These phase bound- 
aries are shown by thick blue lines in Fig. [T] and have a 
ground state degeneracy Wn- Let us concentrate on the 
boundary that separates P2 and P3. The allowed nearest 
neighbor configurations are (+3/2, +3/2), (+3/2, +1/2) 
and (+1/2, +3/2), while the (+1/2, +1/2) is not allowed. 
In the one dimensional chain this rule gives a degeneracy 
Wn — Fff-i + Fn+i, where P/v is the iV-th Fibonacci 
number (W 2 = 3, W 4 = 7, W 6 = 18, W 8 = 47, and so 
on)i2£ In the case of square lattice, we cannot give an 
explicit formula for Wn, numerically we find W% =31 
for the 8-site cluster with P/4 symmetry and W\q = 68 
for the 10-site cluster with C4 symmetry (the degeneracy 
depends on the shape of the cluster). 

Starting from this phase diagram, we study the effect 
of the off-diagonal exchange J below, using perturbation 
theory 



B. Mapping to an effective XX Z model 

Sufficiently far from the A = and h — points, where 
we are essentially dealing with two types of spins only 
(| and I f)), the F3-P2-F1 phase transitions can be 
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FIG. 1: (color online) Phase diagram in the Ising limit as 
the function of the anisotropy and magnetic field. The spin 
configurations on A and B sublattice are shown, as well as 
the degeneracies of the ground state manifolds on the phase 
boundaries (the dashed line is a first order phase boundary). 
Long arrows represent the S* = ±3/2 spin states, while the 
sort ones the S z = ±l/2's. The coordination number £ = 2 
for the chain and £ = 4 for the square. F 1 and F3 are uniform 
phases, while the others break the translational invariance and 
are two-fold degenerate. 



mapped to an effective spin-1/2 model XXZ model: 
n cS = J £ (of of + o\o) + Aof of) -h^t ( 2 ) 



effective interaction terms, and the phases Al and Fl 
correspond to the Ising and the saturated phases of the 
effective model, respectively. 

The effective XXZ model is in a gapped Ising phase 
for A > I . Thus it becomes clear from our mapping that 
the phase P2 disappears once J > J z /3 and the phase 
Al when J > J z /4, with the phase Fl surviving. 

The XXZ-model has been extensively studied in the 
literature, and numerical methods find no trace of super- 
solids on bipartite lattices. Instead, the zero magnetiza- 
tion gapped phase of the XXZ model (P2 in the map- 
ping) is separated by a first order transition from the 
gapless superfluid phase j 37 i 38 The phase separation can 
be prevented, e.g., by longer range diagonal exchanges 13 . 
Likewise, the supersolid phase can also be stabilized by 
introducing second neighbor correlated hoppings (in the 
language of the equivalent hard-core boson problem), 
where the hopping on the second neighbor depends on 
the occupancy of the site along the hopping path i 14 i 25 
Such terms may arise in higher orders of perturbations 
theory, but even then the existence of the supersolids is a 
question of very delicate balance between different terms. 

The physics of the transitions between P2 and PI, and 
PI and ^43 cannot be mapped to an XXZ model in sim- 
ple terms. In that case we shall distinguish sites that can 
be occupied with spins in three different states. Since 
one of the states ("ft) occupies one of the sublattices, and 
the two other states share the the other sublattice, the 
mechanism (see, e.g., Ref. (23[) that leads to phase sepa- 
ration is suppressed and the formation of the supersolid 
is much more natural. 



where the of are the spin-1/2 operators on site i that 
act on the Hilbert space made of the |f) and \\) effec- 
tive spins. Selecting the mapping | ft),] f) — > |f),|l) 
and comparing the matrix elements between the S=3/2 
Hamiltonian ([T]) and the effective Hamiltonian ([2]), we 
obtain the following parameters for the mapping: 



J 
h 



3J' 
3J, 

h-2A- C,J Z 



(3) 

(4) 
(5) 



The Mapping is valid in leading order of the off-diagonal 
exchange. In this case, the P2 phase corresponds to the 
Ising phase of the effective model, and the F3 and Fl 
phases to the saturated phases of effective Hamiltonian. 
Analogously, the mapping | f), | I) — > ||) leads to 



A = ^ 

J 
h 



■h 
47' 
4J, 

h. 



(6) 

(7) 
(8) 



C. Estimating the first order phase transitions 

From the Ising phase diagram we learned that the 
boundary between Al and ^43 is of first order, corre- 
sponding to level-crossing in the energy spectrum that is 
otherwise gapped. We may assume that for not too big 
values of J this holds as well, so that we can estimate 
the corrections to the phase boundary by comparing the 
ground state energies that is expanded in powers of J. 
The lowest order corrections appear in the second order: 



Fax 



A_(J Z 2(J 2 
4 8 



9CJ 2 



N 4 8 (C-l)Ja 32A-8(C + 1)J; 
E A3 = 9A %J Z %J 2 

4 8 (24C- 8)J 2 - 32A' 



(10) 



Comparing these energies, we get that the first order 
phase transition between Al and A2> in the square lattice 
happens when 



4 7 2 

A = 2J Z - + 0(J 4 ) 



(11) 
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for small J. In the case of the one-dimensional chain we 
get 

9 T 2 

A = J z -^-+0(J 4 ). (12) 

Similarly, from the second order corrections given in 
the Appendix, Eqs. (|A4[) and (|A3I) . the boundary be- 
tween the phases PI and Fl is 

A = 2J 2 -^ + 0(J 4 ), (13) 
for a square lattice and 

3 I 2 

A = J z -^- + 0(J 4 ), (14) 

•Jz 

for a chain. 



D. Field induced instability of uniform phases 

The field induced instability of Ising phases can be 
thought of as a softening of magnetic excitations. The 
simplest magnetic excitations corresponds to lowering or 
raising the spins on a site that becomes delocalized due 
to the off-diagonal J term. These excitations are gapped 
in the Ising (plateau) phases, and the value of the gap 
changes with magnetic field and interaction parameters. 
When the energy gap vanishes, it means that these ex- 
citations can be created in arbitrary number and an off- 
diagonal long-range order develops. For small values of J 
we can use perturbation expansion to get the dispersion 
of these excitations. 

In the case of a uniform order the spins on the two 
sublattices are equal, and the perturbational expansion 
of the excitation energy is simple. Let us pick an ex- 
ample, e.g. the instability of the fully polarized phase 
F3 towards the plateau P2. In F3 the ground state is 
Y[j | itj)- A spin excitation in this case corresponds to 
lowering the spin to a f on a given site, with a diagonal 
energy cost 

3 

AE = h-2A- -(J z . (15) 

The off-diagonal terms move the excitations onto the 
neighboring sites, as shown in Fig. [3Ja) , with a 

3 T 

(Mi \n\ Mi) - T (16) 

hopping amplitude, leading to the following dispersion: 
u k = h -2A + ^C(Jlk ~ Jz)- (17) 

Here 

7k = i^exp(ik^), (18) 
^ s 



TABLE II: (color online) Summary of instabilities of uniform 
phases. 
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where the summation is over the vectors 8 pointing to- 
ward the £ nearest neighbors. The quantity 7k takes its 
minimal value —1 at k = (n, . . . ), and its maximal value 
1 at k = (0, . . . ). For the one-dimensional chain (( = 2) 

-f k = cosk x , (19) 

and 

7k = ^ (cos k x + cos k y ) (20) 

for the square lattice (£ = 4). In the F3 phase this ex- 
citation is gapped with a minimum at k = (tt, . . . ), and 
lowering the magnetic field the gap closes when 

h sat = ^((J z + J) + 2A. (21) 

Instabilities of this kind are summarized in Eqs. (|A7jl - 
(|A9[) . the corresponding critical fields are shown in Ta- 
ble HI1 and are plotted in Fig. [^a) for J/J z = 0.2. We 
shall mention that these results are not independent from 
the mapping we discussed in the previous subsection. 

We note that in the case of the F3 phase Eqs. (JT7J) and 
(f2Tj) are exact, while for Fl higher order terms in J/J z 
appear in the dispersion. 



E. Dispersion of spin— excitations in translational 
symmetry breaking states on the square lattice 

The instability (softening) of the excitations in the 
two-sublattice gapped phases (Al, A3, PI, and F2) that 
break the translational symmetry occur in the second 
order of exchange coupling J. Namely, the on-site exci- 
tations on the two sublattices have different energy, and 
depending on the energy difference we shall apply a differ- 
ent scheme for the degenerate perturbation calculation. 
As an example, we discuss the lower instability of the 
2/3-plateau P2 phase. 

The wave function in the Ising limit of the P2 phase 
is given by 

i* p2 >=n n \ (22) 

j&A j'GB 

Applying the S~ operator on the A and the B sublattice, 
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FIG. 2: (color online) (a) Instabilities (thick lines) of the 
gapped phases in the square lattice as obtained from the per- 
turbation theory for J/J z = 0.2. We also show the J = 
phase boundaries of Fig. [T] with thin lines for comparison, (b) 
Variational phase diagram as function of A/£J and h z /(J for 
large exchange anisotropy J/J z = 0.2 and a bipartite lattice 
with £ neighbors. Dashed lines stand for first order, while 
solid lines denote second order phase transitions. The solid 
dot ending the first order transition represent a tricritical 
point. 



respectively. The two energies are identical when A = 
2 J z — this is actually the phase boundary between the 
Pi and F\ phase in the Ising phase diagram, Fig. [TJ 

First we discuss the case when the energy difference is 
larger than J: when AE B - AE A = AJ Z - 2A > J, the 
ground state manifold is given by the \§f) states. Since 
{<&f\H\<&f) = V3J for neighboring i and i' sites and 
{^f\%\^) = 0, the | 4) excitation acquires dispersion 
in a second order process in J, where the | f) excitations 
on the B sublattice can be viewed as virtual state [see 
Fig. 0(b)] . This leads to 



wp 2 ^pi(k) =h-6J z - 



(2) 

where the wp2->Pi denotes additional second order con- 
tributions in J that are independent of k — the full 
form of the dispersion is given in Eq. (| A15|) . In 
other words, the gap closes quadratically for small val- 
ues of J. A similar calculation can be done for the 
AEa — AE B = 2A — 4J 2 3> J case, when the ground 
state manifold is given by the \§f ) states, and we simi- 
larly get a dispersion, Eq. (|A14I) in the Appendix, where 
the hopping amplitude is quadratic in J (we note that 
an additional virtual state assists the hopping). 

When the two excitation have equal energy at A = 2 J z , 
the perturbation theory outlined above obviously fails 
[the hopping amplitudes in both Eqs. (|A14[) and (|Al5j) 
diverge]. In that case we shall include both \$>f) and 
states into the ground state manifold. Actually, 
we can do it also when the energies are not equal, and 
to get the dispersion of the spin excitations, we need to 
diagonalize the following 2x2 problem in k space: 



3J 2 



4 J z - 2A 



.(2) 

J P2^P1 



(27) 



w 



P2 



( h-6J z 4\/3J7k 
y4-y/3J7 k h~2J z ~2A 



(28) 



where we neglected second order contributions. The 2x2 
matrix is easily diagonalized, leading to the 



w k = /i-4J 2 -A± a/(A-2J 2 ) 



48J 2 7 2 



(29) 



dispersion. We notice that for A = 2J Z the dispersion 
becomes linear in J, while for J <C |A — J z \ expanding 
the square root we get 



A-2J Z 



we get 



= Ui>n n itj>ito'>, 

j<i A res 

= iti>n n iti>itfi'}, 



jeA j'<=b 



with diagonal excitation energies 

AE A = h-6J z , 
AE B = h-2A-2J z , 



(23) 
(24) 



(25) 
(26) 



uj k = h - 4 J z - A ± (A - 2 J z ) ± 



(30) 



In other words, we obtain the result of the second or- 
der perturbation theory, Eq. (|27[) . To be consistent, we 
shall take into account all the second order processes that 
contribute to the dispersion. This can be done systemat- 
ically, and the full expression is given in Eq. (|A19[) . The 
critical field at which the gap vanishes can then be de- 
termined without difficulty, and the instabilities of this 
type, given by Eqs. (IA17|) . (|A18I) . and (|A19|) are shown 
in Fig. (Ha) for J/J z = 0.2. 
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FIG. 3: (color online) (a) Schematic figure for the first or- 
der hopping process that occurs during the instability of uni- 
form phases Fl and _F3, where the dispersion is oc 47k. (b) 
Schematic representation of the second neighbor correlated 
hopping that gives the dipersion oc 167^. There are 8 neigh- 
boring places where the magnon can hop through a virtual 
state on the B-site. 



Before proceeding to discussion of the phase diagrams, 
let us mention briefly that for the gapped phases the 
variational wave function is of the same form as it is 
in the Ising-limit, when we neglected the off-diagonal 
terms. Similarly, the expressions for the ground state 
energy are also identical, since to get a contribution from 
the off-diagonal Sf SJ + S?Sj term, we need to tilt the 
spins out of the z axis. Furthermore, the boundary of the 
gapped phases, assuming a continuous phase transition, 
can be determined by studying the stability of the gapped 
variational wave function |^o): the eigenvalue of the 
d 2 E/du a dup indicates a second order phases transition, 
where u a and up are coefficients of a wave functions that 
is orthogonal to l^o)- 



A. Phase diagram in zero magnetic field 



III. VARIATIONAL PHASE DIAGRAM 

In this section we construct the phase diagram vari- 
ationally, assuming either uniform or two-sublattice or- 
dering. We search for the ground state in the following 
site-factorized variational form: 



i*> = n n \^a)Mb) 



(31) 



where 



\i> A ) ex uo| 1r) + e^il t) + ^ 2 u 2 \ |) + e^u 3 \ ty) (32) 

and a similar expression for \tJjb}- In the general case, 
there are 6 independent variational parameters for \iPa) 
and another 6 for |V>b) that are determined by minimiz- 
ing the ground state energy 



E 



(33) 



Recalling that the Hamiltonian is 0(2) symmetric and 
commutes with the S z = J2i &i operator, the state 
rotated by <p around the z axis and given by the 
exp(— ipS z )\^) wave function has the same energy as the 
state described by |\&). We can therefore reduce the num- 
ber of independent parameters for sites A from 6 to 5, so 
in total we have reduced the number of independent pa- 
rameters from 12 to 11. It appears, however, that all the 
phases we have found are coplanar, and after a suitable 
rotation the amplitudes in the wave function can all be 
chosen to be real. 

The site-factorized variational wave function (1311) is 
actually indifferent to the connectivity of the lattice, the 
only information about the lattice that enters the expres- 
sion of the energy is the number of the neighbours £. For 
concreteness, we look at the case of the square lattice, 
however the results can be easily generalized to any bi- 
partite lattice by replacing J — > C,J/A and J z —> (Jz/4 in 
equations and phase diagrams shown below. 



First let us take a look at the zero field phase diagram. 
We find two - the completely and the partially aligned 
- axial antiferromagnetic states, A3 and Al as well as 
a superfluid U(l) phase between them. This latter is 
referred to as a planar state in Ref. [39[ for the spins 
are aligned in the lattice plane, but also can be called 
superfluid since spin-rotation symmetry breaking phases 
exhibit finite spin stiffness^! that is the property of such 
phases. In the following we refer to this phase as SFo. 
Between the planar superfluid SF and the two axial an- 
tiferromagnets Al and A3 an additional superfluid phase 
appears. The in-plane components of this conical antifer- 
romagnet have the same properties as the planar super- 
fluid but it exhibits finite staggered magnetization too, 
inheriting the property of the antiferromagnetic phases. 
Therefore we call this phase SFa- A schematic figure of 
the various phases is shown in the phase diagram in Fig. 

E " 

The relevant order parameters for zero external field 
are the staggered magnetization m s z = 
the superfluid order parameter Ojjn 
where S+ = (Sf,S^). O 

staggered magnetization. The expectation values of or- 
der parameters as a function of A/J z are shown in Fig. 
[5] for various values of J/J z . 

The first order phase boundary between the two axially 
aligned antiferromagnetic phases is A = 2J Z . It can be 
determined by comparing the ground state energies of A3 
and Al listed in Table HI 

The ground state wave functions of sites A and B in 
the planar superfluid phase SFq can be expressed as 
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Sg\, and 



- ilS-L 



,. ~f7(i) is actually the in-plane 



= e^ s i|* SF ), 

with a single variational parameter rj: 
1 
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FIG. 4: (color online) Variational phase diagram for ft = as 
the function of A/J z and J/J z . Solid lines stand for contin- 
uous (second order) phase boundaries, while the dashed line 
denotes the first order phase boundary of the phase A3 . 
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FIG. 5: (color online) Order parameters as the function of 
A/J z for different values of J/J z - The axial antiferromagnets 
have finite staggered magnetization, and zero U(l) order pa- 
rameter. In the planar superfluid phases mj* = but the 
expectation value of Otr(i) is finite, while the conical antifer- 
romagnet exhibits both type of order. 



ip can take arbitrary value and is related to the U(l) 
symmetry breaking (we recall that the Hamiltonian com- 
mutes with the S z operator), as it determines the direc- 



tion that the spins point to in the xy plane 

377(77 + 1) 



(*a\S a \* a ) 



3?7 2 
377(77 



3t/ 2 



cos ip, 
sm<p, 



0. 



(37) 

(38) 
(39) 



The ground state energy as the function of parameter 77 
reads 

E$ F ° (V) _ 3 V 2 + 3_ A _ I877 2 (v + lf j (4Q) 



N 



4 3?7 2 



1 



(St? 2 



1 



In the energy expression the J z term is absent, as this 
wave function has only spin component in the xy plane. 
Minimizing the energy gives a cubic equation for 77. How- 
ever, a given 77 value minimizes the energy for the 



A _ 3(r/ 2 - 1)(3t7+ 1) 
J ~ 



3?7 2 



parameter in the Hamiltonian. 
find 



1 



(41) 



For small values of A we 
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A 2 



0(A 3 ) 



144J 2 

and the ground state energy can be approximated as 



(42) 



n 



3 A 2 

-A-— + 0(A 3 ) 

4 16J 1 ' 



(43) 



that gives the phase boundary with the antiferromagnetic 
phase A3 



J=Jz 



A 
3" 



A 2 
72J Z 



0(A 3 ) 



(44) 



as seen in Fig. 2J 

For A = 0, when the anisotropy is absent, 77 = 1 and 
Eq. dSH]) is just a spin coherent state of the spin-3/2 
Neel-state of the SU(2) symmetric Heisenberg model ro- 
tated into the xy plane. For A > the S z — ±3/2 
components in the wave function are suppressed. In the 
A — > +00 limit the 77 = A/3 + 0(1) and we are left 
with a wave function with S z = ±1/2 spin components 
only. Out of these two states we can mix a spin point- 
ing to arbitrary direction, however the length of the spin 
is not constant - it is the largest when lying in the xy 
plane (the length is then 1 as opposed to 1/2 when point- 
ing in the z direction, a consequence that they are still 
5 = 3/2 spin). For this reason the antiferromagnetic ex- 
change term gains the most energy with the planar spins, 
as in Eq. (|3l>|) . When the exchange interaction becomes 
anisotropic, and the SfS z term becomes strong, this en- 
ergy can compensate the directional length dependence 
of the spin, and can choose a spin configuration with a 
finite z and xy component. This happens in the conical 
superfluid phase, denoted by SF A in Fig. @] 
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The phase boundary of the conical superfluid phase 
(SFa) towards the planar phase (SF ) and fully polar- 
ized AFM phase (A3) is a complicated expression. It is 
shown in Fig. 2J Starting from phase Al at a given A 
value, a second order phase transition occurs to the su- 
perfluid phase SFa- When the exchange coupling J is 
large enough, in-plane spin components appear continu- 
ously as we reach into SFa- The ground state can be 
expressed as it follows: 

\* A ) oc e- ll ^(|^)+ W |t)+"|l)+H^)) (45) 
|* fl ) cx e- ,;(v+7r) ^(u;|^+«|t>+"ll> + l4»(46) 

with real u, v and w variational parameters. 

The instability of the partially aligned AFM phase Al 
against canting gives the phase boundary 



J 



J Z (J,-A) 
J 2 -4A 



(47) 



between Al and SFa- 

The same model for one dimension has been treated 
by mean field calculations in Ref. [39| for quantum spin 
1/2, 1 and 3/2. The phase diagram for the case 5 = 3/2 
is similar to our findings, however the conical superfluid 
phase SFa is missing due to a more restricted variational 
wave function they used. 



B. Heisenberg exchange with on— site anisotropy 

In the following we discuss the phase diagram as the 
function of magnetic field and single-ion anisotropy when 
the exchange between two neighboring sites is SU(2) sym- 
metric (i.e. the J = J z Heisenberg model with on-site 
anisotropy) . The phase diagram outlined in Fig. [5] was 
calculated by the variational method introduced above. 

On the h z = line the ground state is the planar super- 
fluid phase (SFq) introduced previously. As the magnetic 
field becomes finite, the spins cant out of the plane con- 
tinuously, and the superfluid ground state SFf exhibits 
finite magnetization m z alongside the finite staggered in- 
plane order parameter Of/(i) [Fig. [7]- A schematic figure 
of the conical SFf is shown in Fig. [6] and the ground 
state wave function can be characterized as: 

|*A> « e-'^i (| H) + u\ t) +v\i)+ w\ 11)) (48) 
|* B ) cx e-^+^Bd^+uW+vlD+wlH-))^) 

where u, v, and w are all real numbers. 

The ground state energies of the axial ferromagnetic 
phases and the fully polarized ferromagnetic state are 
given in Table HI 

Analytical expression for the ground state energy of 
SFf is beyond our reach, however, the phase boundary 
with the neighboring _F3 phase can be given by calcu- 
lating the critical field for the polarized phase. This is 
exactly the same as the instability approximation for _F3 




FIG. 6: (color online) Phase diagram in the function of A/ J 
and h z /J (Jz = J). 



a/c;j z =o 




F3 



Ou ( i; 



1/3 
1 

1/2 




A/y z =3/2 






F3 



1/3 

1 


A/g 2 =3 


SF F 




1/2 


















1/3 

E? 1 
1/2 




a/c;j z =9/2 
-_SF F 


F1 


sf f ; 









4 6 



10 



FIG. 7: (color online) Order parameters as the function of 
magnetic filed for different values of A parameter. The fully 
and partially polarized ferromagnets F3 and Fl exhibit finite 
magnetization m z , while in the conical ferromagnetic phase 
SFf the expectation values of both m z and Odti) are finite. 



in the case of the Ising limit, and is given by the same 
Eq. (|21l) . Above the saturation field the fully polarized 
ferromagnetic phase is stabilized. For large enough values 
of A the spins become shorter and the partially polarized 
plateau phase Fl emerges. Calculating the instability of 
Fl, the phase boundary turns out be 



h = J + 2 J z + A ± a/ J 2 - 14JA + A 2 . 



(50) 
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As expected, from the mapping to the effective XXZ 
model (Sec. IIIBI) . we found no evidence for gapped 
phases that break the translational symmetry. 



C. The effect of exchange anisotropy and the 
emergence of supersolid phase. 

Finally let us examine the collective effect of exchange 
and single-ion anisotropics as well as the magnetic field. 
In the previous subsection we learned that only the fer- 
roaligned spins in PI and P3 are present as gapped 
phases for the case of the Heisenberg exchange (J z = J), 
with a superfiuid phase (canted antifcrromagnet) in be- 
tween them. As the value of J/J z is lowered, islands 
of plateaus and antiferromagnetic phases emerge in the 
sea of the superfiuid phase. We choose a relatively large 
anisotropy J/J z = 0.2, as in that case we learned from 
the perturbational expressions that the 2-fold degenerate 
gapped phases might be stable, as shown in Fig. (Ha). In- 
deed, the variational phase diagram, shown in Fig. HJb) , 
displays all the phases we were looking for: The super- 
fluid phase takes place around the axial ferromagnets, 
while between the plateaus and axial antiferromagnetic 
phases - i.e. the gapped phases that exhibit staggered 
diagonal magnetic order - a very robust supersolid phase 
arises. 

The extension of the supersolid around the phases Al 
and P2 is the broadest at their tips, when A is not too 
large. As we increase A, the supersolid region decreases, 
and eventually vanishes for A — > +00. Since in this limit 
the mapping to the XXZ model becomes exact fSec. HIBj) . 
our finding is also consistent with numerical works on the 
XXZ model on the square lattice that do not seem to find 
supersolid i 13 ' 37 i 38 

The variational calculation finds all the phase bound- 
aries to be of second order, except a single first order one 
around A « 2 J z [shown in dashed line in Fig. [5{b)] that 
is inherited from the J — phase diagram, Fig. [TJ 

The expression of the phase boundaries of the axial 
ferromagnetic phases are the same as in the Heisenberg 
limit (see Eqs. (f2"Tj) and (|5T))0 . We determined the phase 
boundaries of the plateaus and axial antiferromagnetic 
states by calculating spin wave instability. We found that 
the boundary for the 2/3 plateau can be given as 



A = - - 6 ± v/(9-3J z ) 2 -9J 2 , 



(51) 



with 6 = 2 J z + h /2-3j ■ Similar calculations give 



h 
£ 



y/2{Jl - J 2 ) + 2{J Z - A) 2 - 2S, 

y/(J 2 - A(A - 2J Z )) 2 + 32 J 2 A(2A - J z ), (52) 



for the phase boundary of the partially polarized axial 
antiferromagnetic phase Al, 
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FIG. 8: (color online) Expectation value of order parameters 
per site as the function of h/J z for different values of A/J z . In 
the axial antiferromagnetic phases Al and A3 only the stag- 
gered magnetization has finite expectation value. For A = 
there is a first order phase transition from the completely 
polarized A3 phase to the superfiuid phase. All the other 
field induced transitions are second order transitions. The 
superfiuid phase exhibits finite magnetization m z and finite 
staggered in-plane magnetization Ojy(i). The ferromagneti- 
cally ordered phases F3 and -Fl are characterized by finite 
magnetization m z , while the plateaus (PI and P2) have an 
additional finite staggered magnetization m s z . In the super- 
solid phase all four order parameters have finite expectation 
values. 



with 8 = 2J Z + 
PI plateau, and 



6, J z 
h/2-3J z 



for the phase boundary of 



A = 3 J z - 



4 



9J 2 



(54) 



A = 0- -±-y/(0 — 3J Z 



9J 2 



(53) 



for the boundary of the axial antiferromagnetic phase 
A3. When J = Eq. (J53j) and (JM) give back the h = 
6J Z — 2A phase boundary that separates A3 and PI in 
the Ising limit. The ground state energies and phase 
boundaries for the superfiuid and supersolid phases can 
only be obtained numerically. The ground state wave 
function for the superfiuid with ferromagnetic m z is given 
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by Eq. (|49)) , and for the supersolid by 

\* A ) oc e-^(\il)+u\t)+v\i)+uU)) (55) 
\* B ) oc e -^+-)^(|^)+ u '|t)+^U)+^'|4)X56) 

where u, u', v, v' , w, and it)' are all real. Fig. [5] shows 
the evolution of the order parameters which can be used 
to find out the nature of the phases as we increase the 
magnetic field for a few selected values of A/J z . 



D. S=l phase diagram 

At this stage, it is useful to compare the predictions 
of the variational approach to a better studied problem. 
Supersolid phases have been found in spin-1 anisotropic 
Heisenberg antiferromagnet in Ref. l23l . so we constructed 
the variational phase diagram for this model as well. The 
Hamiltonian is identical to Eq. [TJ but now with S = 1 
spin operators. The phase diagrams, for vanishing J 
and J = 0.2J Z , are shown in Fig. [HI In the Ising limit, 
Fig. [HJa), we find two uniform phases (denoted as 00 
and 11, using the values of the S z components) and two 
phases breaking the translational symmetry: the 11 with 
zero magnetization and the 10 one-half magnetization 
plateau. The XXZ-like physics can be identified for the 
transition between the 11,10, and 00 phases, where the 
supersolidity is a fragile phase. The region between the 
10 and 11 is of different nature, and we expect the su- 
persolid to be robust in this part of the phase diagram. 
And that is exactly the region where Ref. found su- 
persolidity. Furthermore, the nature of the phase tran- 
sitions is also in qualitative agreement, inasmuch the or- 
der of the phase transitions is concerned. Specifically, 
we recover the first order transition between the upper 
boundary of the 10 phase and the superfluid. It is also 
useful to compare Fig. HJb) to the phase diagram of the 
one-dimensional chain obtained by DMRGi 2 ^ the extent 
of the gapped phases is reduced in the chain, and the 
supersolid survived only in a small region close to the ll 
phase. 

The calculation of the phase diagram is quite straight- 
forward - assuming two sublattice order, the variational 
wave function is given by Eq. pip , now with 



3/2 



1/2 



11 


(a) 


\ 10 








11 ' 


00 


col 





1/2 



11 




(b) 


/ SF > 








10 






xSS y 




11 


1 / 






1 
1 
1 
1 


00 



FIG. 9: (color online) The phase diagram of the anisotropic 
S=l model in the (a) Ising-limit for a bipartite lattice with 
coordination number C, and (b) for the square lattice ((" = 4) 
when J = 0.2J Z , obtained from the variational calculation. 



phase boundaries 

h 2 = ((J z - A — (J) (C J z - A + C J), (58) 
h 2 = A(A-2CJ), (59) 
[h - A)(CJ 2 + A - h)(h - QJ Z + A) = 2C 2 J 2 A,(60) 

respectively. 



IV. SUPERSOLID IN THE ONE-DIMENSIONAL 
MODEL 



oc + e^uolO) + e^mll) (57) 

and a similar expression for \iPb}- We are now dealing 
with 8 independent variational parameters altogether, 
that can be reduced to 7 by using the U(l) symmetry 
of the model. Similarly to the spin-3/2 case, we get so- 
lutions where all the amplitudes can be chosen to be real 
numbers for the Hamiltonian we look at. 

The saturation field is given by h sat — A + QJ Z + (J, 
and from the stability analysis of the 11, 00, and 10 
gapped phases we get the following equations for their 



In this section we complement the variational study 
using a variant of the Density Matrix Renormalization 
Group^ (DMRG) method on the anisotropic 5 = 3/2 
spin chain. A Quantum Monte-Carlo study has shown 
that a supersolid phase can realized in the anisotropic 
5 = 1 spin chain j 2 ^ a result confirmed by DMRG calcu- 
lations in Refs. I28l - l30l Therefore it is plausible that a 
supersolid states is also present in the anisotropic spin- 
3/2 chain. 

We map out the phase diagram for the chain and search 
for signatures of supersolid phases. The DMRG method 
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we used optimizes variationally a wavefunction based on 
a matrix-product stated (MPS) Ansatz for an infinite 
chain. Algorithms using this approach are efficient in 
one dimensional systems because they exploit the fact 
that the ground-state wave functions are only slightly 
entangled. For mapping out the phase diagram, we used 
comparably small matrix dimensions of \ = 50, while 
for estimating the central charge we used matrices up to 
X = 200. 



m s, z entanglement entropy 
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FIG. 10: (color online) Zero-field phase diagram for the infi- 
nite chain as a function of A/J z and J/J z , as obtained from 
DMRG calculation with x = 25. Panel (a) shows the stag- 
gered magnetization. Panel (b) shows the half chain entangle- 
ment entropy, i.e., the von-Neumann entropy of the reduced 
density matrix for a bipartition of the chain into two half 
chains. 

The zero magnetic field phase diagram is shown in 
Fig. EH We can clearly identify the gapped A3 and Al 
uniaxial phases with finite value of the staggered magne- 
tization to|' and small entanglement entropy, and a gap- 
less phase with algebraic correlations (Luttinger liquid). 
The extension of the Al phase is limited to J/J z < 0.25 
values, following the estimate based on the mapping to 
the XXZ model in Sec. HTbI 

Fig. [IT] shows the phase diagram in the present finite 
magnetic field. Again, the gapped phases can be identi- 
fied using the uniform and staggered magnetization (m z 
and m S z), and the small entanglement entropy. The ex- 
tension of the gapped phases essentially follows the vari- 
ational phase diagram (see Fig. HJb)). However, the su- 
persolid phase is more fragile in the one-dimensional case 
due to strong quantum fluctuations. Consequently, the 
gapless phase in the phase diagram is predominantly a 
simple Luttinger liquid with algebraically decaying cor- 
relations and characterized by the integer central charge 
that measures the number of gapless modes. We calcu- 
lated the central charge of the gapless phases using the 
method outlined in Ref. 42] for a few selected points in 
the phase diagram. Within numerical accuracy we find 
c = 1, as shown in Fig. 1121 This is in accordance with 
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FIG. 11: (color online) Phase diagram as the function of A/J z 
and h/J z for (a-c) J/.J Z = 0.1 and (d-f) J/J z = 0.2. We 
show the uniform and the staggered magnetization along the 
z axes, where the plateau phases can be identified. The large 
increase of the entanglement entropy indicates gapless phases. 
The phase diagram is obtained from DMRG calculation with 
X = 25. 

our expectation originating from the mapping to the ef- 
fective XXZ model, that the gapless phases between the 
F3 and P2, P2 and Fl, and Fl and Al are all Luttinger 
liquids. 
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FIG. 12: (color online) Estimate of the central charge from 
the entanglement entropy for four different points in the LL 
phase and one point in the supersolid phase (the h/J z = 1.95, 
A/J z — 0.5 point). The solid lines are fits based based on the 
5=1 In 5 + const, formula, with c set to 1. In all these points 
the gapless phases are characterized by c = 1 central charge. 

We have searched for a supersolid phase in the vicin- 
ity of the gapped phases that break the translational in- 
variance. We made a scan by varying the field for a 
fixed value of A/J z and J/J z ; the results are plotted in 
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FIG. 13: (color online) The magnetic field dependence of 
the order parameters as a function of the magnetic field for 
J/J z = 0.2 and A/J z = 0.5. The curves are result of DMRG 
calculations with x = 40 and with a small field h x /J z — 10~ 4 . 
In (a), the non-vaninshing off-diagonal order parameter m% 
shows the extension of the gapless phases. The finite value of 
the m a z and m x indicate a robust supersolid phases, as seen 
in (b). 



Fig. Q21 For the simulations, we added a tiny magnetic 
field of order 10 -4 along the x-axis to break the U(l) 
symmetry around the-z axis. A finite value of the di- 
agonal (staggered magnetization mf) and off-diagonal 
(magnetization along the x axis, m x ) order parameter 
indicates the presence of the supersolid. It appears that 
the supersolid is stable in a small region only, between 
the A3 and PI gapped phase, with a continuous phase 
transitions. Both the magnetization and the staggered 
magnetization in the supersolid show a square root like 
behavior at the lower and upper critical fields, like the 
magnetization in XXZ model does, see for example in 
Ref. [43j . This is due to the density of the states of the 
spinons, and is already observed in the XY model, when 
it is mapped to free fermions. Recall that the density of 
states of free fermions has a van Hove singularity at the 
band edges, and this shows up as a square root singular- 
ity in the magnetization curve of XY (and XXZ) model. 
In Fig. [14] we straighten out this singularity. This singu- 
larity is also inherited for the staggered magnetization at 
the critical fields. The central charge in the supersolid 
is also c = 1 (the lowest line in Fig. [T^J here we set the 
h x = 0, as otherwise a finite h x induces a gap in the 
spectrum) . 

From variational calculations, we expect a continuos 
phase transitions into the supersolid also at the upper 
edge of the P\ phase. Numerically, however, we find a 
first order transition into the LL phase. 




FIG. 14: (color online) The magnetization has a square root 
singularity at (a) lower h c ,i/J z = 1.8579, and (b) upper 
h c ,2/Jz = 2.0195 critical field. The solid lines show the 
m\ « 2.68(ft - h e ,\)/Jt + I6.9(h - h c ,i) 2 /J% and (1 - m z f w 
1.27(/i c ,a — h)/J z +13.2(h c ,2 — h) 2 / J z fits to the magnetization 
curves. 



V. EXACT DIAGONALIZATION STUDIES 

To get further insight into the problem in higher di- 
mensions, we have numerically diagonalized small (8- and 
10-site) clusters of spin 5 = 3/2 arranged on the square 
lattice with periodic boundary condition and searched 
for signature of different phases in the energy spectrum. 
The two-sublattice states break translation symmetry, 
so we expect two degenerate ground states with momen- 
tum k = (0,0) and (7T, tt) in the thermodynamic limit. 
In the gapped phases, these two levels are well separated 
from the other states, while in the supersolid, where both 
translational symmetry and U(l) symmetry breaking oc- 
curs, we expect two copies of the Anderson towers in 
the spectrum that is the signature of the U(\) symmetry 
breakingi 44 i 45 Unfortunately, the large spin makes the fi- 
nite size scaling difficult, and without a finite size scaling 
we cannot be sure about the exact nature of the ground 
state. Nevertheless, even our small cluster gives impor- 
tant support for the variational phase diagram. 

In Fig. [T5] we show the energy spectrum for the C4 
symmetric 10 site cluster and J = 0.2 J z around A = 2J Z , 
where we expect the first order transition from the phase 
A3 into the supersolid to happen. In zero field the 
ground state has S z = 0, and in Fig. [l~5Tb) we see that 
the energy level curvatures of lowest lying states in the 
k = (0,0) and (tt,tt) sector are essentially indistinguish- 
able for A < 1.88 J z and well separated from the higher 
levels. This indicates the presence of a gapped, two- 
sublattice state that we can associate with the A3 phase. 
The sharp level anti-crossing at A w 1.88J Z indicates 
a first order transition. In the S z — 1 sector we ob- 
serve the spin excitations with a narrow bandwidth and 
a k «-» (tt, tt) — k symmetry, following Eq. (|A16I) as cal- 
culated from the perturbation theory in Sec. Ill El For 
A > 1.88J Z , the k = (0,0) and (tt, tt) levels are also 
close, and the these two levels are equally close and re- 
versed in order for S z = 1 in Fig. I15f a). an indication for 
the U(l) symmetry breaking, possibly with translational 
symmetry breaking (the supersolid phase). 
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FIG. 15: (color online) The first few lowest lying energy levels 
of a 10 site cluster for (a) S z — and (b) S z = 1 as a function 
of A/J z . We set J = 0.2J Z . The inset shows the available k- 
points in the Brillouin zone. 
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FIG. 17: (color online) Magnetization as a function of mag- 
netic field, as obtained from variational calculation and exact 
diagonalisation. Here J = 0.2J Z , and A/J z =0, 1.5, and 4.5. 
ftsat is the saturation field [Eq. (I21[) ]. 
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FIG. 16: (color online) The gap A = E^^) — i?(o,o) as a 
function of A/J z and h/J z for the 10-site cluster. The solid 
curves separate the different S z sectors. 



In Fig. [16] the energy gap between the k = (0, 0) and 
(7r, 7r) ground states in the different S z sectors is shown as 
a function of A/ J z and magnetic field, as this may serve 
as an indicator of the translational symmetry breaking. 
We can identify the gapped phases (except for Al) and 
their extension is even quantitatively in good agreement 
with the variational phase diagram, shown in Fig. Hfb). 
The consistency between the variational and exact di- 
agonalization result is also supported in Fig. 1171 where 
we compare the magnetization calculated by these two 



VI. CONCLUSIONS 

In this paper we studied the effect of exchange and 
easy-plane anisotropics on the formation of magnetiza- 
tion plateaus and supersolids in spin-3/2 system on (un- 
frustrated) bipartite lattices, with the aim to extend the 
results of earlier studies on the stability of supersolids 
in anisotropic spin-1 model on the square lattice^ and 
spin-1/2 bilayer system a 22 ' 24 ' 25 to larger values of spins. 

In the Ising limit (J = 0) we find both uniform and 
translational symmetry breaking magnetic phases with 
gapped excitation spectrum with zero of finite mag- 
netization (magnetization plateaus). We discussed the 
macroscopic degeneracy of the ground state at the phase 
boundaries and showed that when the off-diagonal ex- 
change interaction J becomes finite this degeneracy is 
lifted and new gapless phases emerge. All the plateaus 
continuously evolve from the Ising limit, and the degen- 
eracy of the boundaries in the Ising-limit gives a hint on 
the order of the phase transition and on the nature of 
the gapless state. Not surprisingly, our variational calcu- 
lation shows that the supersolid phases are concentrated 
around the plateaus that break the translational symme- 
try. In particular, the tendency toward supersolidity is 
greatly enhanced when the degeneracy of the boundary 
is 2 x 2 N / 2 (due to the choice of two states on the sites 
of one of the sublattices, while the sites of the other sub- 
lattice are occupied with a third type of states), as in 
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this case the diagonal translational order is preformed, 
and the off-diagonal order is easily established on the 
sublattice occupied with the two states. In addition, for 
large anisotropies we have confirmed the stability of the 
plateau states using perturbation theory, and found a 
good agreement between the two approaches regarding 
the extension of the gapped phases. 

In zero field we plotted the variational phase diagram 
as function of the on-site and exchange anisotropies. 
Aside from the axial antiferromagnetic phases and planar 
supcrfluid phase, we find a biconical superfluid which si- 
multaneously exhibits the diagonal and off-diagonal stag- 
gered characteristics of the former phases. 

In the J = J z Heisenberg limit, when the exchange 
interaction is SU(2) invariant (but we keep the on-site 
anisotropy A that breaks the SU(2) symmetry), the 
plateau and supersolid phases disappear and only the 
uniform phases and the superfluid phase between them 
are present. 

The variational phase diagrams for zero and finite mag- 
netic field were compared to DMRG calculations carried 
out in one dimension. We have found convincing evidence 
for a supersolid state that is realized in a region between 
two gapped phases that break the translational symme- 
try. For a two-dimenal square lattice, we performed ex- 
act diagonalisation in two dimension on small clusters for 
J/J z = 0.2 and identified the characteristics of various 
phases from the energy spectrum. The extension of the 
gapped phases based on these calculations proved to be 
in good qualitative and quantitative agreement with the 
variational findings. 

Our study was initially inspired by the Ba2CoGe2C>7 
layered material, where Ref. [35| estimates A/J z w 8 and 
J w J z . While these anisotropies are not strong enough 
to stabilize a magnetization plateau, an anomaly occurs 
around m/m sa t = 1/3 in the magnetization curve - this 
is also observed experimentally: the magnetization curve 
changes it's slope at h « 9T4& 
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Appendix A: Perturbation expansion 

Here we are presenting the results of the Rayleigh- 
Scrodinger perturbation theory applied to states and ex- 
citations in the J — > limit. 

1. Second order corrections in J to the 
ground-state energy 

The second order correction to the energy/ (per site) of 
the different phases are as follows: 



,(2) 
~-Al 



.(2) 
'-A3 



.(2) 
'-F1 



.(2) 
'PI 



.(2) 
'PI 

-(2) 
'F3 



8J 2 _ 9J 2 
3J Z 2(4A - 5 J z ) 

9J 2 
2(11J Z -4A) 
12J 2 



2A - J z 

6J 2 
7 J z - 2A 

2 J z 



(Al) 

(A2) 

(A3) 

(A4) 

(A5) 
(A6) 



2. First order degenerate perturbation theory for 
excitation spectrum of the uniform Fl and F2 phases 



ujfi^P2 = -h + 2 J z + 2A + 6 J7 k 
Wfi-s-AI = h-2J z + 8Jju 
^F3^P2 = h — 6 J z — 2A + 6 J7k 



(A7) 
(A8) 
(A9) 



where 7k is defined in Eq. 



3. Second order degenerate perturbation theory 
for excitation spectrum of the staggered phases 
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oc t2 q t2 aq t2 

— = " + 2A ~ eJ > - 5^55 - 3iray 16 ^ + T^ffi (A10 » 

UJ w . + g + ? «L c _ E ^ ic(1 » )S + 8 ) (AH) 

27 J 2 12 J 2 64J 2 36 J 2 2 J 2 2 , 3 J 2 a , 

— - = ^ + 2J -4A3IX-2A^X + ^T + iA^5X-— ( 16 ^ + 8 )-2A3IX 16 ^ (M2) 

q/2 /2 

wp 2 -vf 3 = -h + 6J z + 2A- — (l6 7 2 + 8) + 12— (A13) 

5J 2 J z 

— 1 = 71 " 2J * " 2A " 2J7T2A - 8X 16 ^ - 2A3IX 16 ^ + X (AW) 

21 j2 3 j2 

— - = /l - 6 ^ + IJT " iJ-^A 167 ^ (A15) 

1 O r2 oc 72 q t2 

— ' " -" + 6J ' - 2A - ToJT^A + ilf^A - sIST^' 16 ^ + 8) (A16) 

In case that we include the excitations on both sublattices, the S~ excitations from the Al in the k space are 
eigenvalues of the 

^ " I 473J 7k 2A - 2 J, - A - - .^(le. 2 +8) J" ^ 

matrix. If we expand for J up to second order, this will give Eq. IA121 the corrections to the dispersion directly to the 
Fl phase. 

Similarly, for the PI phase 



! 6J 7k 2J z + / l -2A-|^----i M (16 7 2 +8) 1 



and for the P2 phase: 



-6 J 2 + h - 4V3J 7k | _ ^(2) 

4^3 J 7k -2 J 2 + h - 2A - - f£ (lfryg + 8) 
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